Boundary conditions for Einstein's field equations: Analytical and numerical analysis 
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Outer boundary conditions for strongly and symmetric hyperbolic formulations of 3D Einstein's 
field equations with a live gauge condition are discussed. The boundary conditions have the property 
that they ensure constraint propagation and control in a sense made precise in this article the 
physical degrees of freedom at the boundary. We use Fourier-Laplace transformation techniques to 
find necessary conditions for the well posedness of the resulting initial-boundary value problem and 
integrate the resulting three-dimensional nonlinear equations using a finite-differencing code. We 
obtain a set of constraint-preserving boundary conditions which pass a robust numerical stability 
test. We explicitly compare these new boundary conditions to standard, maximally dissipative 
ones through Brill wave evolutions. Our numerical results explicitly show that in the latter case 
the constraint variables, describing the violation of the constraints, do not converge to zero when 
resolution is increased while for the new boundary conditions, the constraint variables do decrease 
as resolution is increased. As an application, we inject pulses of "gravitational radiation" through 
the boundaries of an initially flat spacetime domain, with enough amplitude to generate strong 
fields and induce large curvature scalars, showing that our boundary conditions are robust enough 
to handle nonlinear dynamics. 

We expect our boundary conditions to be useful for improving the accuracy and stability of 
current binary black hole and binary neutron star simulations, for a successful implementation 
of characteristic or perturbative matching techniques, and other applications. We also discuss 
limitations of our approach and possible future directions. 



I. INTRODUCTION 



In many numerical simulations in General Relativity one integrates Einstein's field equations on a spatially compact 
domain with artificial timelike boundaries, effectively truncating the computational domain. This raises the question 
of how to specify boundary conditions. In this article we address this question within the context of the Cauchy 
formulation, in which the field equations split into evolution and constraint equations. We adopt a free evolution 
approach, in which the constraints are solved on the initial time slice only and the future of the initial slice is computed 
by integrating the evolution equations. Boundary conditions within this approach should ideally satisfy the following 
three requirements: (i) be compatible with the constraints in the sense of guaranteeing that initial data which solves 
the constraint equations yields constraint-satisfying solutions, (ii) permit to control, in some sense, the gravitational 
degrees of freedom at the boundary, and (iii) be stable in the sense of yielding a well posed initial-boundary value 
problem (IB VP). 

There are several motivations for the construction of boundary conditions satisfying the above properties. First, 
recent detailed analysis of binary neutron star evolutions 1] showed that the presence of artificial boundaries in 
current state of the art evolutions, which use rather ad hoc boundary conditions, dramatically affects the dynamics 
in the strong field region, near the stars. While it can be argued that this effect should disappear when placing 
the boundaries further and further away from the strong field region, until ideally, the region of interest in the 
computational domain is causally disconnected from the boundaries, in practice, this would require huge computer 
resources, especially in the three-dimensional case. Even though some kind of mesh refinement should help in placing 
the boundaries far away, there is in any case a minimum resolution needed in the far region in order to reasonably 
represent wave propagation, thus constraining the size of the computational domain for a given amount of memory. 
Next, an understanding of boundary conditions is important if the evolution equations include elliptic equations, since 
in this case effects from the boundaries can propagate with infinite speed, and so have an immediate effect on the 
fields being evolved. Examples of cases in which elliptic equations arise include elliptic gauge conditions (such as 
maximal slicing @| or minimal distortion conditions 'sf) and constraint projection methods. Finally, isolating 

the physical degrees of freedom at the boundaries should be important in view of Cauchy-characteristic (see for 
a review) or Cauchy-perturbative ^ Il0| matching techniques, where a Cauchy code is coupled to a characteristic 
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or perturbative "outer module" which is well adapted to carry out the evolution in the far zone. In these cases, it 
is important to communicate only the physical degrees of freedom at the boundary between the two codes, since the 
Cauchy code and the outer module might be based on completely different formulations. 

Boundary conditions satisfying requirement (i) can be constructed by analyzing the constraint propagation system, 
which constitutes an evolution system for the constraint variables and is a consequence of Bianchi's identities and 
the evolution equations. If it can be cast into first order symmetric hyperbolic form, the imposition of homogeneous 
maximally dissipative boundary conditions [TH Il2| for the constraint propagation system guarantees that a smooth 
enough solution of the evolution system (if it exists) which satisfies the constraints initially automatically satisfies 
the constraints at later times. Homogeneous maximally dissipative boundary conditions consist in a linear relation 
between the in- and outgoing characteristic fields of the system, which is chosen such that an energy estimate can be 
derived. This energy estimate implies that the unique solution of the constraint propagation system with zero initial 
data is zero, i.e. that the constraints are preserved during evolution. Maximally dissipative boundary conditions 
for the constraint variables usually translate into differential conditions for the fields satisfying the main evolution 
equations, since the constraint variables depend on derivatives of the main variables. This means that the resulting 
boundary conditions for the main system usually are not of maximally dissipative type and, as discussed below, 
analyzing well posedness of the corresponding IB VP becomes more difficult. 

Requirement (ii), controlling the physical degrees of freedom, is a difficult one since there are no known local 
expressions for the energy or the energy flux density in General Relativity. Nevertheless, one should be able to control 
the physical degrees of freedom in some approximate sense, as for example in the weak field regime approximation, in 
which one linearizes the equations around flat spacetime. In this approximation, it might be a good idea to specify 
conditions through the Weyl tensor, since it is invariant with respect to infinitesimal coordinate transformations of 
Minkowski spacetime. More precisely, since there are two gravitational degrees of freedom, we should specify two 
linearly independent combinations of the components of the Weyl tensor. Below we will discuss boundary conditions 
which involve the Newman-Penrose complex scalars and \E'4 (see, for instance, with respect to a null tetrad 
adapted to the time-evolution vector field and the normal to the boundary. If the boundary is at null infinity these 
scalars represent the in- and outgoing gravitational radiation, respectively. Furthermore, it turns out that these scalars 
are invariant with respect to infinitesimal coordinate transformations and tetrad rotations for linear fluctuations of any 
Petrov type D solution represented in an adapted background tetrad 0| . This class of solutions not only comprises 
flat spacetime but also the family of Kerr solutions describing stationary rotating black holes. Since in many physical 
situations one is interested in modeling asymptotically flat spacetimes, such that if the outer boundary is placed 
sufficiently far away from the strong field region spacetime can be described by a perturbed Kerr black hole, we 
expect the boundary condition \l/o = to be a good approximation for a "non-reflecting" wave condition. These 
boundary conditions are actually part of the family of conditions imposed in the formulation of Ref. [Tsl l , to date the 
only known well posed initial-boundary value formulation of the vacuum Einstein equations, and were also considered 
in [3. 

Requirement (iii), the well posedness of the resulting IB VP, turns out to be a difficult problem as well: For 
quasilinear symmetric hyperbolic systems with maximally dissipative boundary conditions there are well-known well 
posedness theorems J_8j J^J which state that a (local in time) solution exists in some appropriate Hilbert space, 
that the solution is unique, and that it depends continuously on the initial and boundary data. The proof of Ref. 
|15| is based on these techniques. There, using a formulation based on tetrad fields, the authors manage to obtain 
a symmetric hyperbolic system by adding suitable combinations of the constraints to the evolution equations in 
such a way that the constraints propagate tangentially to the boundary. In this way, the issue of preserving the 
constraints becomes, in some sense, trivial. (See |2(]| for a treatment in spherical symmetry in which the constraints 
propagate tangentially to the boundary as well.) However, for the more commonly used metric formulations it seems 
difficult to achieve tangential propagation for the constraints with a symmetric or strongly hyperbolic system, and 
therefore, one has to deal with either constraint propagation across the boundary or systems that are not strongly 
or symmetric hyperbolic. Here we choose to deal with constraint propagation across the boundary and strongly or 
symmetric hyperbolic systems. There has been a lot of effort in imdcrstanding such systems, both at the analytical 
UL22.M,2£2h,M2L2&.l^, 30] and numerical ^ 16 , 21, 21 M, M, M level. Although partial proofs of weU 
posedness have been obtained using symmetric hyperbolic systems with maximally dissipative boundary conditions 
|22. 123I I24L I23 , it seems that these kind of boundary conditions are not flexible enough since constraint-preserving 
boundary conditions usually yield differential conditions. 

In this article we construct constraint-preserving boundary conditions (CPBC) for a family of first order strongly 
and symmetric hyperbolic evolution systems for Einstein's equations. For definiteness, we focus on the formulation 
presented in Ref. |34t| , which is a generalization of the Einstein-Christoffel type formulations |3^ |3^ IstI Ib^ with a 
Bona-Masso type of gauge condition 39, 40] for the lapse. However, our approach is quite general and should also 
be applicable to other hyperbolic formulations of Einstein's equations. In Section ^1 we briefly review the family of 
formulations considered here and recall under which conditions the main evolution equations are strongly or symmetric 
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hyperbolic. In Section ITTll we discuss the constraint propagation system and analyze under what circumstances it is 
symmetric hyperbolic. Having cast this system in symmetric hyperbolic form we impose CPBC via homogeneous 
maximally dissipative boundary conditions for the constraint variables. These conditions are differential boundary 
conditions when expressed in terms of the fields satisfying the main system. We complete these boundary conditions 
in Section llVl bv imposing extra conditions which control the physical and gauge degrees of freedom. There are several 
possibilities for doing so, of which we analyze the following two: 1) we first specify algebraic conditions in the form 
of a coupling between the outgoing and ingoing characteristic fields (referred to as CPBC without Weyl control later 
in this article); 2) we specify boundary conditions via the Weyl scalars ^I'o and 5*4, as discussed above (referred to as 
CPBC with Weyl control later in this article) . 

Section [V] is devoted to an analysis of the well posedness of the resulting IBVP. As mentioned above, this is a 
difficult problem since our boundary conditions are not in algebraic form. I particular, they are not in maximally dis- 
sipative form, so we cannot apply the standard theorems for symmetric systems with maximally dissipative boundary 
conditions. Therefore, our goal here is more modest: We analyze the IBVP in the "high frequency limit" by con- 
sidering high-frequency perturbations of smooth solutions. In this regime the equations become linear with constant 
coefficients, and the domain can be taken to be a half plane. In this case solutions can be constructed explicitly by 
performing a Laplace transformation in time and a Fourier transformation in the spatial directions that are tangential 
to the boundary. This leads to the verification that a certain determinant is non-zero. If this determinant condition 
is violated, the system admits ill posed modes growing exponentially in time with an arbitrarily small growth time. 
Therefore, the determinant condition yields necessary conditions for well posedness of the resulting IBVP and allows 
us to discard several cases which would lead to an ill posed formulation. These ill posed formulations appear even in 
cases in which both the main and the constraint propagation system are symmetric hyperbolic. We also stress that the 
determinant condition verified in this article is a weaker version of the celebrated Kreiss condition 0, li^ , which 
yields well posed IBVP for hyperbolic problems with algebraic boundary conditions |4ll l44l. However, the Kreiss 
condition guarantees nothing for the present case of differential boundary conditions. See |30j for an example of an 
IBVP with differential boundary conditions which satisfies the Kreiss condition but fails to be well posed in L^. 

In section IVII we discretize the IBVP by the method of lines. The spatial derivatives are discretized using finite 
difference operators that satisfy the summation by parts property and we explain in detail how we numerically 
implement our CPBC. Then, in sections IVIII and IVIIII we perform the following 3D numerical tests: First, we evolve 
IBVPs which fail to fulfill the determinant condition and are therefore ill posed. The goal of evolving these systems 
is to confirm the expected lack of numerical stability. We do confirm such instability, finding an obvious lack of 
convergence: the results exhibit exponentially in time growing modes, where the exponential factor gets larger as 
resolution is increased, as predicted by our analytical analysis. Next, we focus on evolutions of two systems that 
do satisfy the determinant condition, differing only on whether they control components of the Weyl tensor at the 
boundary or not. Finally, we compare stable evolutions of systems with CPBC with evolutions where maximally 
dissipative boundary conditions are given for the main evolution system, and are therefore expected to violate the 
constraints. 

We first analyze evolutions of these four systems through a robust stability test |45ll46| |. In this test, random initial 
and boundary data is specified at different resolutions, and the growth rate in the time evolved fields is observed. A 
growth rate that becomes larger as resolution is increased is a strong indication of a numerical instability, while for 
numerical stability the growth rate should be bounded by a constant that is independent of resolution. We find that 
systems that violate the determinant condition fail to pass the robust stability test, as expected. However, we also 
find that at least some systems with CPBC with Weyl control which satisfy the determinant condition are numerically 
unstable as well, although in a somehow weaker sense (explained in the text) that reminds the numerical evolution of 
weakly hyperbolic systems j^. In contrast to this, the systems with CPBC without Weyl control which satisfy the 
determinant condition that we have evolved pass the robust stability test for the length of our simulations (usually 
between 100 and 1,000 crossing times). 

Next, we concentrate on evolutions of Brill waves |48j . and confirm the expectations drawn from the robust stability 
test in what concerns numerical stability. Using these waves we further concentrate on a detailed comparison between 
the results using maximally dissipative boundary conditions for the main evolution system and stable CPBC. Our 
convergence tests strongly suggest that in the former case the constraint variables do not converge to zero in the limit 
of infinite resolution, implying that one does not obtain a solution to Einstein's field equations, while in the latter 
case for the same resolutions the constraint variables do converge to zero. 

Next we concentrate on the stable CPBC case and evolve pure gauge solutions, using high order accurate finite 
difference operators which satisfy the summation by parts property. The operators used are eighth order accurate in 
the interior points and fourth order accurate at and near the boundary points. 

As a final numerical experiment, we also concentrate on the stable CPBC case and inject pulses of gravitational 
radiation through the boundaries of an initially flat spacetime. We inject pulses of large enough amplitude to create 
very large curvature in the interior (as measured by curvature invariants), showing that our CPBC are strong enough 
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to handle very non-linear dynamics. The order of magnitude achieved by the curvature invariant measured corresponds 
to being at roughly r — 0.7 from the singularity, in a Schwarzschild spacetime of mass one. In these simulations this 
curvature is produced solely by the injected pulses. 

A summary of the results and conclusions are presented in Sect. IIXI Technical details, like the derivation of the 
constraint propagation system and of the characteristic fields, and a special family of solutions to the linearized IB VP 
with Weyl control arc found in 1X1 IbI and O 



II. STRONGLY HYPERBOLIC FORMULATIONS WITH A LIVE GAUGE 

In this section we review the family of hyperbolic formulations of Einstein's field equations constructed in |34l |. 
which is an extension of the Einstein-Christoffel type of formulations which incorporates a generalization of 

the Bona-Masso js^^] slicing conditions. It consists of 34 evolution equations for the variables {N, gij, Kij,dkij , Ai}, 
where N is the lapse function, gij the three-metric, Kij the extrinsic curvature, and where the extra variables dkij 
and Ai represent the first order spatial derivatives of the three-metric and of the logarithm of the lapse, respectively. 

The evolution equations are obtained from the ADM evolution equations in vacuum by adding suitable constraints 
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to the right-hand side of the equations (see El for more details). Following the notation of Ref. |3J|, we have 

doN = -F{N,K,xn, (1) 
dog^, = ~2K,^ , (2) 
doK,, = R,, - + T%Ak - A,A, - 2K^ 'Ki, + KK,, 

+ 19vC + Cg'''Cki^J)l , (3) 

2 

dodkrj = -2dkK,j - 2AkKij +rigk(iCj) +X9ijCk + jj9i{idj)dkl3'' , (4) 
doA. = -^A.~-—{g-%Kk,-d.uK-^)--—^+^C.. (5) 

Here, do = {dt — £p)/N, and F{N, K,x^) is a function that is smooth in all its arguments and that satisfies a — 
(OkF) / (2N) > 0. For the simulations below we shall choose F — N ■ K which corresponds to time-harmonic slicing, 
but for our analytical results we shall leave F unspecified for generality. We assume that the shift vector /3* is a fixed, 
a priori specified vector field. The parameters |6^ 7, ^, 77, x, ^ control the dynamics off the constraint hypersurface, 
defined by the vanishing of the following expressions 

C = l{g^iRt,i-K'^iKki+K^), (6) 

a = {dkKu - d,Kki) + ^(d' - 2h'')Kk^ + \d^'''Ku , (7) 

Ckij = dkij — dkgij , (8) 
Cikij = d[idk]tj , (9) 
C|^' = A,~d,N/N, (10) 



^If ^ (11) 

Here, C = is the Hamiltonian constraint, Ci = the momentum one, and Ckij = 0, Cikij = 0, =0 and 

^if^ ~ ^ artificial constraints that arise as a consequence of the introduction of the extra variables dkij , Ai |63j | . 
The Ricci tensor Rij belonging to the three-metric is written as 

^ij — 2^*°' {—dkdiij + dkd(^ij-ji + d(^^d\ki\]) — dad^^ki) 

+ ^U,ki + \{du~ 2bk)T% ~ T%T\^ , (12) 
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where hj = dkijg'^^ dt = dkijg^-' and 
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The evolution equations H1I2I3I4I5|I have the form of a quasilinear first order system, 

dou = A{uydjU + T{u), 



(13) 



where u — (N, gij , Kij , dkij , Ai) and the matrix- valued functions A-', j — 1,2,3, and the vector- valued function J- 
are smooth. We are looking for solutions with given initial data on some three-dimensional manifold which satisfies 
the constraints equations. In order to guarantee the existence of solutions we restrict the freedom in choosing the 
parameters a, 7, C, "Hi X-, £, by demanding that the corresponding initial-value formulation is well posed. This means 
that given smooth initial data, a (local in time) solution should exist in some appropriate Hilbert space, be unique, 
and depend continuously on the initial data. Although in this article we are interested in the numerical evolution of 
spacetimes on a spatially compact region with boundaries, well posedness of the problem in the absence of boundaries 
is a necessary condition for obtaining a numerical stable and consistent evolution inside the domain of dependence of 
the initial slice. An easy and intuitive way of finding necessary conditions for well posedness is to look at high frequency 
perturbations of smooth solutions: Let p = {t,x^) be a fixed point and Uq a smooth solution in a neighborhood of p. 
Perturb uq according to u{t,x^) = UQ(t,x^) + eu{ujt) exp{iu!njX^), with e, uj real, and rii a constant one-form on R'^ 
which is normalized such that g{pY^ninj = 1. If we evaluate the evolution equations at a point near p, divide by toe, 
and take first the limit uj ^ 00 and then the limit e ^ 0, the evolution equations reduce to 



dtu^ [N{p)A{p,n)+(3{pynj]u, 



(14) 



where the matrix A{p^ n) is given by 



A{p, n)u 



—2aniK - 



K) 



- X 9i] {Kkr. 

n,K) 



UkK) 



(15) 



where we have set A-,- = bj — dj and where the index n refers to the contraction with n'. Notice that by rescaling 
and rotating the coordinates a;* one can always achieve that g{p)ij — Sij, and by rescaling the coordinate t one can 



achieve that N{p) = 1. For this reason, we drop the entry p in the following. We call the system H14|l the associated 
frozen coefficient problem. A necessary condition for the well posedness of the initial-value formulation defined by 
the (non-linear) Eqs. (|1I2I3I4I5|) is the well posedness of the associated frozen coefficient problem. If some extra 
smoothness properties are satisfied (see|Hll this condition is also a sufficient one. The frozen coefficient problem is well 
posed if the matrix NA{n) + P^rii is diagonalizable and has only real eigenvalues for each n. Clearly, this is true if and 
only if the matrix A{n) is diagonalizable and has only real eigenvalues. As we have shown in :34j this can be easily 
analyzed by taking advantage of the block structure of A{n): Suppose u is an eigenvector of A{n) with eigenvalue fi. 
Then 

,(0) 

where — {N,gij), m^^^ = (Kij), m'^^ — {dkij,Ai), and where the matrices A and B are read off from Eq. (|15(l . A 
sufficient condition for A{n) to be diagonalizable and posses only real eigenvalues can be obtained by considering the 
equation 




Explicitly, we have 



/x^ify = ABK^j 



y?Kij = Kij + An(^^Kj)n + Bmn-jK + Cgij {Knn - K) , 



(16) 



where the coefficients A, B and C are 



A = -2-|- J(3C-l)r7-C, 
B = | + i(3C-l)r; + e+l + 2a, 



(17) 
(18) 
(19) 
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We now demand that AB is diagonahzable and has only positive eigenvalues. As we have shown in [3J| this guarantees 
that A{n) is diagonahzable and has only real eigenvalues. Representing AB in an orthonormal basis e], ef, such 
that ej = rii, we find 

tJ.^K„„ = il + A + B)Kr,n + {B- C)Ki , (20) 

fi^KnA = {l+A/2)KaA, (21) 

^^^Ki - (1 - 2C)Ki , (22) 
fi^KAB = Kab , (23) 

where A — 2,3 and Kab — Kab — SabS'~''^ Kcd/'^- From this we immediately see that AB is diagonahzable with 
only positive eigenvalues if and only if 

Ai = 2a, (24) 
A2 = i(27 + l)(2 + 2x-ry)-ic^, (25) 

A3 = -Ic^x-v + ^O'lcv, (26) 

are positive and if B — C = Ai + ^(A2 + 1) — 2A3 = whenever Ai = A2. In^we derive the characteristic fields, which 
are given by the projections of u onto the eigenspaces of A{n). These fields play an important role in the construction 
of boundary conditions. Using these fields, we also derive in^sufficient conditions for the non-linear evolution system 
to be strongly hyperbolic and thus yield a well posed initial value formulation. 



III. CONSTRAINT PROPAGATION SYSTEM 



In order to obtain a solution to Einstein's equations, not only do we have to solve the evolution equations but 
also the constraints. We want to follow a free evolution scheme, in which the constraints are solved initially only. 
For this scheme to be valid, we have to guarantee that any solution for such initial data automatically satisfies the 
constraints in the computational domain everywhere and at every time. At the numerical level, since the constraints 
are already violated initially due to truncation or roundoff errors, we have to guarantee that the numerical solution to 
the evolution equations converge to a constraint-satisfying solution to the continuum equations in the limit of infinite 
resolution. In order to show this, the constraint propagation system, which gives the change of the constraint variables 
under the flux of the main evolution system, plays an important role. In this section we show that for a suitable range 
of the parameters 7, Xi ^ this system can be cast into first order symmetric hyperbolic form. We then specify 
boundary conditions that guarantee that zero initial data for this system leads to zero constraint variables at later 
times. 

The constraint propagation system is derived in^ up to lower order terms which are linear algebraic expressions 
in the constraint variables and whose precise form are not needed for the purpose of this article. In order to analyze 
under which conditions the system is symmetric hyperbolic, it is convenient to decompose Cikij into its trace and 
trace-less parts, 

Cikij = Eikij + - {gi(iB.j)k - gk(tB.j)i) + - gijWik , 

where Eikij = E^ik](^ij) is trace-less with respect to all pair of indices and where Bij and Wij are determined by the 
traces Sh = C^^i^^, Ah = C\^^]^ and Vik = Ci^^/, 

4 12 4 9 12 

Bki = 77 Ski — ^ Aki — - Vki , Wik = - Vik + — Aik ■ 
3 5 5 5 5 

(Notice that Sij is symmetric trace-less while Aij and Vij are antisymmetric.) In terms of the variables U = 

{Cl^\C'kij, Eikij, C,Ci, SijAij,Vi.j,clf^), where Ay = Ay + F„72 + Cfj^\ the non- linear constraint propagation 
system has the form 



doU = Ac{uydjU + B{u)U, 



(27) 
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where the principal symbol Ac{u,n) = Ac{uynj is given by 



Ac{u, n)U 



( 



-i(2 + 2x-??)C„ 
-(27 + l)niC + C'S'm-i„ 

[i(2x-?7)+C] "[.C,] 
(j] + 3x) n[,Cj] 



J 



and where the matrix B{u) depends on the main fields u and their spatial derivatives, but not on U. The system Eq. 
(|27|l is called symmetric hyperbolic if there exists a symmetric positive definite matrix H which may depend on u but 
not on n such that llAciu,n) is symmetric for all u and n. From the above representation of the principal symbol 
it is not difficult to see that the system is symmetric hyperbolic if the following conditions hold 



(27 + 1) (2 + 2x - ?7) > 0, 

Cv<o, 

2x - ?/ + 4^ < 0. 



(28) 
(29) 
(30) 



Notice that these conditions automatically imply that A2 > and A3 > 0, which are necessary conditions for the main 
evolution system to be strongly hyperbolic. If the conditions H28I29I30() are satisfied, a symmctrizer is given by the 
quadratic form 



2 + 2x - ?7 3 



Z A^i A ■ ■ 



(31) 



where 



A, 



A^, +-Vi 



a 



(A) 



A. 



a 



{A) 



2-V ■ 

4(ry + 3x) 
(2x - ??) + 4^ 

(2x - ??) + 4^ 



A,. 



The symmetrizer allows us to obtain an energy estimate for solutions to Eq. (|27|) on a domain O of K'^ with smooth 
boundary dO and suitable boundary conditions on dO. Defining the energy norm 



E. 



constraint 



ts^ / U^UUd^x, 



(32) 



differentiating with respect to t and using the constraint propagation system, Eq. (|27|l . we obtain 



at 



2 / L/^H [N{Aljd,U + BU)+ P^diU] 
Jo 

I U'^nNAc{n)U (fx 

JdO 

[ [NUB + NB'^U - d^{NB.A'c + H/3'')] U d^x, 
Jo 



(33) 



where here rii denotes the unit outward one-form to the boundary. In the last step we have used the fact that Ac (n) is 
symmetric with respect to the scalar product defined by H31|l and assumed that the shift is tangential to the boundary 
at do. As one can easily verify. 



U^JIAc{n)U = ^5'-' ( 



J 1' 3 



(34) 
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where T^/^'' = —Ci ± {Ani — C,Sni) ± (27 + Ijn^C are the in (+) and out (— ) going characteristic constraint fields. 
Therefore, if we impose the boundary conditions 

= , (35) 

where the matrix [S^) satisfies g^-' S^viVk < g^^Vj^vi for aU one-forms v^, we obtain the energy estimate 

~^-^constraints — COTiSt • E^Q^g^j- dints : {^^) 

where the constant only depends on bounds for NB and H^-^diiNHAlj + H/3*). Since Econstraints > this proves 
that Econstraints{t) = for alH > provided that Econstraints{i = 0) = 0. For this reason, we call the three conditions 
constraint-preserving boundary conditions. 



IV. BOUNDARY CONDITIONS FOR THE MAIN EVOLUTION SYSTEM 

In this section we consider the main evolution system, Eqs. H1I2I3I4I5|I . on a open domain O of R'^ with smooth 
boundary dO. We also assume that the shift vector is chosen such that it is tangential to dO at the boundary. This 

means that at the boundary we have six ingoing characteristic fields, denoted by , "^ab' ''^aai '^ii? (see|Blfor their 
definition; here and in the following. A, B = 2,3, and quantities with a hat denote trace-free two by two matrices), 
and thus we have to provide six independent boundary conditions. Following the classification scheme of Ref. one 
can show that the first field is a gauge field, the second ones are physical fields, and the last are constraint-violating 
fields. We stress that this classification scheme does only make precise sense in the linearized regime for plane waves 
propagating in the normal direction to the boundary (see |3 for a more detailed discussion about this). 

If we forgot about the constraints, we could give data to the six ingoing fields. The simplest possibility would 
be to freeze the ingoing fields to their values given by the initial data. Provided the evolution system is symmetric 
hyperbolic this would yield a well posed IBVP. However, in the presence of constraints, the boundary conditions have 
to ensure that no constraint-violating modes enter the boundary, i.e. the boundary conditions have to be compatible 
with the constraints. In fact, the numerical results of section rVIIII show explicitly that freezing of the ingoing fields 
to their initial values does not, in general, provide a solution to Einstein's equations: The constraints are violated. 

Instead of the freezing non-constraint preserving boundary conditions just mentioned, we impose the three 
constraint-preserving boundary conditions (|35|) . This fixes three of the six conditions we are allowed to specify 
at do. We complete these three conditions in the following two ways: 

1. CPBC without Weyl control 

Here we adopt the simplest possibility and impose algebraic boundary conditions on the "gauge" and "physical" 
fields, 

v^n^J = aviJ+h, (37) 

^AB = , (38) 

where a and b are constants satisfying \a\ < 1, |6| < 1, and h and Hab are functions on dO describing the 
boundary data. The justification for the bounds on a and b will become clear in the next section. Notice that 
the choice a = b ~ —1 results in Dirichlet conditions in the sense that data is imposed on some components of 
the extrinsic curvature while the choice a — b ~ 1 yields boundary conditions that impose data on combinations 
of spatial derivatives of the three metric (see the definitions of vim and 

v'^AB inlH). In our simulations below, 

we choose a = b = which yields Sommerfeld-like conditions. 

2. CPBC with Weyl control 

Here we replace Eq. (|38ll by a similar condition for the Weyl tensor. We impose the boundary conditions 

viV = avU+s, (39) 
w'^AB = cw^ab+Hab , (40) 

where |c| < 1 and where w^ab defined in terms of the electric {Eij) and magnetic (i?y) parts of 
the Weyl tensor in the following way: Let ei, 62, 63 be an orthonormal triad with respect to the three 
metric gij such that ei coincides with the unit outward normal to the boundary. Then, w^ab ~ 
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6^6^/2 j (^Eij ± n'^ekiiBj'-) , where Skij is the volume element associated to gij. If the vac- 
uum equations hold, the electric and magnetic components can be determined by 

E^j = doK.j + K.aK'^j + V(,A,) + AAj , (41) 
5y = ersi^y'^K',■) , (42) 

where one uses the evolution equation (jSJ in order to reexpress d^Kij in terms of spatial derivatives of the main 
variables. The boundary conditions H4Q(I correspond to the conditions imposed by Friedrich and Nagy 15]. In 
the symmetric hyperbolic system considered in 15j, where the components of the Weyl tensor are evolved as 
independent fields, these boundary conditions arise naturally when analyzing the structure of the equations 
since they give rise to maximally dissipative boundary conditions. In particular, a well posed initial-boundary 
value problem incorporating the condition (|40ll is derived in In contrast to this, the boundary conditions 
l|4U|l are not maximally dissipative for our symmetric hyperbolic system since they are not even algebraic. 

The conditions (|40|l can also be expressed in terms of the Newman-Penrose scalars and with respect to 
a null tetrad which is adapted to the boundary in the following sense: Let eg denote the future-directed unit 
normal to the t = const slices. Together with the above vectors ei, 62, 63 it forms a tetrad. From this, we 
construct the following Newman-Penrose null tetrad: 

l>^ = ^ (e^ -f ei") , k'-^ (eo - 4) , - ^ {4 + ^4) ■ 



Then, we find 



1^22 "'■'"23 ' — "^22 T ""23 ' 

and the boundary condition (|40|l is simply 

*o = c*:, (43) 

where the star denotes complex conjugation (we could generalize this boundary condition by allowing for complex 
values of c). Notice that and ^"4 are not uniquely determined by the unit normals eo and n: with respect 
to a rotation of 62, 63 about the angle 0, these quantities transform through ^ e^^'^^E'o, ^'4 1— > e~^*'^\E'4, 
the factor 2 reflecting the spin of the graviton. However, the boundary condition (|43|) is indeed invariant with 
respect to such transformations. In our simulations below, we shall choose c = corresponding to an outgoing 
radiation condition. 

Finally, as discussed in the introduction, ^'q and 5*4 represent, respectively, in- and outgoing radiation when 
evaluated at null infinity and are gauge-invariant quantities for linearizations about a Kerr background. 



V. FOURIER-LAPLACE ANALYSIS 



We now analyze the well posedness of the IB VP defined by the evolution equations (|1I2I3I4I5|I . the CPBC H35|) 
and the boundary conditions and or and (j^ . We derive necessary conditions for well posedness 
by verifying a certain determinant condition in the high frequency limit. We assume that the parameters are such 
that the evolution equations are strongly hyperbolic since otherwise the problem is ill posed even in the absence of 
boundaries, and such that the constraint propagation system is symmetric hyperbolic. 

Let p € do be a point on the boundary. By taking the high frequency limit we obtain the associated frozen 
coefficient problem at p. After rescaling and rotating the coordinates if necessary, we can achieve that N{j)) = 1, 
g{p)ij — 5ij and that the domain of integration is > 0. In this way, we obtain a constant coefficient problem on 
the half space. Introducing the operator do = dt — f3^di, it is given by 

doK,, = i {-^''dk^J + (1 + O^'^rffe-jfe + (1 - 0^i^bJ) - - 2% A,)) 

+ |%a'=(5fe-4) (44) 

dodkzj = -2dkK,j +T]gk(i{d'^Kj)i - dj^K) +xgi] {d'^Kki - dkK) , (45) 

doA, = -2adJ< + e {d'K,, - d,K) . (46) 
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Notice that this system is equivalent to the one that one would obtain by linearizing the evolution equations around 
flat spacetime in a slicing with respect to which the three metric is flat, the lapse is one and the shift is constant and 
tangential to the boundary, but not necessarily zero. 

Since we have a linear constant coefHcient problem on the half plane, we can solve these equations by means of a 
Laplace transformation in time and a Fourier transformation in the and directions (42„ .43j . That is, we write 



, x^^ — u{ijjx^^ exp(w(2:t + iu:A{x'^ + where 



the solution as a superposition of solutions of the form uit^x^ ^x" 
z G C with Re(z) > 0, S-'^^Coa^b ~ 1, A = 2, 3 and u £ L^(R+). (Notice that for such solutions, dou = ujzu.) 
Substituting this into Eqs. H44I45I46|I one obtains a system of ordinary differential equations coupled to algebraic 
conditions. Since there are six in- and six outgoing modes, there are twelve independent differential equations. The 
remaining equations which are algebraic can be used in order to eliminate the characteristic fields which have zero 
speeds, and one ends up with a closed system of twelve linear ordinary differential equations. Because the system 
is strongly hyperbolic we expect [111 lil] exactly six linearly independent solutions that decay as x^ ^ oo, and six 
solutions that blow up as X"'^ ^ oo. Since we require the solution to lie in we only consider the six decaying 
solutions. The determinant condition consists in verifying that the boundary conditions with homogeneous boundary 
data annihilate these six solutions. If the determinant condition is violated, the problem admits solutions of the form 
u{t,x^,x'^,x^) = u{ujx^) exp{uj{zt + iiUA{x^ + P^t))) for some Re(2:) > 0, where ui can be arbitrarily large, and the 



system is ill posed since \u{t,x^, 



|m(0, x^,: 



'■^)\ — exp{ujKe{z)t) ^ oo as cj ^ cx3 for each fixed t. Thus, 



the determinant condition is necessary for the well posedness of the IBVP and, as we will see, will yield nontrivial 
conditions. 

A convenient way for finding the six decaying solutions is to look at the second order equation for Kij, which is a 
consequence of Eqs. (|44l45l4(j|) . 



(47) 



where the coefficients A, B and C are given in Eqs. (|17I18I19|I . We show below that there are exactly six solutions 
to this second order system which have the form Kij{t,x^ ,x^ ,x^) — Kij^ujx^) exp{Lo{zt + iuiA{x"^ + P^t))), with 
Re(z) > and such that Kij decays as —> oo. Since in Fourier-Laplace space the operator do is just multiplication 
with the nonzero factor loz, corresponding solutions to the original system can be obtained by determining dkij and 
Ai from Eqs. and (|^ . respectively. 

In order to find the decaying solutions of Eq. H47|l we make the ansatz Kij = kij eyipiuovx), where here and in the 
following we set x = x^ for simplicity, and where v \s & complex number with negative real part to be determined. It is 
also convenient to introduce a unit two- vector i)^ which is orthogonal to loai and to decompose kij in the components 



ansatz, Eq. (|17|l splits into the following two decoupled systems 



kABi^^r}^ and k^^ 



rj rj )kAB/'^- Using this 



h 



( hxx ^ 

I.A 

"-A 


= M2{v) 


'^A 

kxLO 


\ ki^uj / 




\ kuiui / 



(48) 



where the matrices Afi(z^) and M2{i') are given by 



Xsv^-l \Av 



fXiiy^ + C-l {B-C)iy^ + ^ i{A + 2C)v 
2C-B ■ " ' " ~ 



X2iy'' ~ ^ - B + C -1 iiA + 4C)iy -{A + 2C) 



V 



,2 _ 4 

^{A + 2B)i^ \{A + AB)v 

B ~\{A + 2B) 



A 
2 



-c 

+ 2 
Av 

1-1/ 



and Ai, A2, A3 are defined in Eqs. (|24I25I26|I . 

The matrix M\ (y) has the eigenvalue-eigenvector pair 



2 2 1 

z = V — \, 



z' = A3(z.2-l), 
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The two vectors are always linearly independent from each other since Re(z) > 0. This yields the solution 



XT] 



LUJ] 



CTi 



1 

— ZI^4 



(72 



1 



(49) 



where cti, a2 are two constants and where here and in the foUowing vi = \J ^z^ + 1, / = 1,2,3,4 (A4 = 1) where 

the branch of the square root for which Ke{i'i) > for Re(z) > is chosen. Similarly, after obtaining the eigenvalues 
and eigenvectors of M2{i') one obtains the solution 



where 



k 



1 



_ —UJUaX , _ —UJUsX 

<J3V3 e + cTiVi e 



aeve e 



(50) 



^^3 



V5 



-1 



\ 



-11^3 

iv3 

2{B - C)vl +A + B + 2C 
2{A + B + 2C)vl +A-B + 4:C 
-i//2(A + 3B) 
-i(A + 3B) 



if Ai ^ A2, U5 = 




if Ai = A2 



Therefore, we have obtained six linearly independent solutions which decay exponentially as x — > cxd and thus lie 
in L^(0). They are parameterized by the constants (which depend on uj and z) ui,...,<jq. A necessary condition 
for the IBVP to be well posed is that these constants are uniquely determined by the boundary data. Before 
checking this condition, it is instructive to have a closer look at the six-parameter family of solutions given by 
Eqs. (|49|) and H5U|I . Let us first compute the momentum constraint variable Ci = Kij ~ diK: It has the form 
Ci(t,a;\x^,a;^) = Ct{ujx^) exp{uj{zt + iwAix^ + P^t))) where 



Cx = -^c^(l-z^|)a4e-"''-'^ + 26^^1/2^5 6-""'^", 



2iLU(f^ e 



(51) 
(52) 

(53) 



where (T5 = (1 — i/|)(Ai — A2)<J5 if Ai 7^ A2 and = (75 otherwise. Thus, Ci = for this family of solutions if and only 
if (72 = 174 = CT5 = 0. In other words, the three-parameter subfamily of solutions parametrized by CT2, (74 and (75 are 
constraint-violating modes. Next, consider an infinitesimal coordinate transformation parametrized by a vector field 
(X^) = (/, X'), and assume zero shift for simplicity. With respect to such a transformation, the linearized lapse and 
extrinsic curvature change according to 



N^N + dtf, 
K,j i-> K,j - d.djf. 

On the other hand, the linearization of Eq.(^) around a Minkowski backgroundl^^ yields 

dtN = -2a K. 



(54) 
(55) 



(56) 



We see that the choice f{t,x^,x'^,x^) = ui ^ exp{uj{zt — vix -f iujAX^)) leaves Eq. (|56|l invariant and induces the 
transformation 







( kxx 




k^ 




k^ 

^A 




k 


1 — > 


k 




\k^^ j 






) 



-1 

—ii/i 

_i 
2 



(57) 
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while Kxri and K^^ri remain invariant. Therefore, it is possible to gauge away the solution parametrized by erg and 
we call this solution a gauge mode from hereon. The remaining family of solutions parametrized by Ci and (T3 are 
physical modes: They satisfy the constraints and cti and are gauge- invariant. 

Next, we verify that the integration constants cti,...,CT6 are uniquely determined by the boundary conditions. First, 
we notice that the expressions (|51I52I53|) for the Fourier-Laplace transformation of the momentum constraint yield a 
three-parameter family of solutions for the constraint propagation system, Eq. H27|l . Since this system is symmetric 
hyperbolic and since we specify homogeneous maximally dissipative boundary conditions for it (see Eq. (|35|l 'l. the 
corresponding IBVP is well posed. In particular, zero is the only solution with trivial initial data. This implies that 
a2 — (Ti — = 0(6^. We stress that such a conclusion cannot be drawn if the constraint propagation is strongly but 
not symmetric hyperbolic, see Ref. p5| for a counterexample. 

Next, using Eqs. H45|) and (|46|l . we find the following expressions for the relevant characteristic fields at the boundary 



,(±) 



-'ab 



1 



Ai ujz 



[XAK + c(i - n)d^KA:c + Ai(i - n)c^] , 



= Kab T 



w[^j^ = lozKab T 



d:,KAB - (1 + Od(AKB) 
TF 

dxKAB - d(AKB)x 



TF 



— [-2cjdAdBK + S,di^ACB)]^^ : 
UJZ 



(58) 

(59) 

(60) 
(61) 



where 



n = 



1 + 2Ai + A2 ~ 4A3 
2(Ai-A2) 



if Ai 7^ A2 and Q arbitrary otherwise. 



(62) 



and where [...]"^^ denotes the trace-free part with respect to the metric Sab- Plugging into this the six-parameter 
family of solutions given in Eqs. H49|l and H5()|l . taking into account the vanishing of the constraint-violating modes, 
cr2 = o'4 = (T5 = 0, and evaluating at x = 0, we obtain 



~(±) 



~(±) 



Z 



Ai 



± 



± z^z^ + 1 



^Vz^ + l ± 



;Vz2 + l ± 



^ +2 



0-3 + 0-6 T 



C 



Ai z 



(1/40-3 + j/io-g) 



03 , 



03 



CTl 



2z 



{z =F J/iC] 0-6 



(63) 
(64) 
(65) 
(66) 
(67) 



In particular, the fields and are gauge- invariant. This is to be expected since these fields are linear 

combinations of the components of the Weyl tensor, which is a gauge-invariant quantity in the linearized regime. The 
are gauge-invariant too, but the remaining fields depend on the constant erg which represents a pure gauge 



fields 



degree of freedom. In particular, if the parameters of the formulation are such that one can choose f2 = 1, are 
pure gauge fields. Now it is not difficult to check the determinant condition for the different boundary conditions 
proposed in the previous section. 



A. CPBC with Weyl control 



Here we analyze the family of boundary conditions H37I4U|I . where we set the boundary data h and Hab to zero. 

We first analyze the boundary conditions w'ab — ^''^ab^ where |c| < 1: Plugging into this the expressions from Eqs. 
(|65I67|I . we obtain 



(Z + Vz2 + 1)2 + C(Z - Vz2 + 1)^ 



(z + + \f - c{z - y/z^ + ry 



01 =0, 
03 =0. 



(68) 
(69) 
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Since the function z i-^ {z + yzM-T)^ maps Re(z) > onto the outside of the unit disk minus the negative real 
axis, while the function z i— > (z — y/z^+l)'^ maps Re(2:) > onto the inside of the unit disk minus the negative real 
axis, and since |c| < 1, it follows that the gauge-invariant constants ui and vanish. Now it also becomes clear why 
we restricted the range of the parameter c: If |c| > 1, there are nontrivial exponential growing solutions with either 
CTi ^ or 0-3 ^ 0, and the system is ill posed. 

In order to analyze the consequences of the boundary condition vit} — a,vix\ we assume that either ft = 1 (in 
which case wit'' are pure gauge fields) or = — Ai. In both cases we obtain, using 0-3 = and Eq. Ht)3|) . 

(76 =0. 



[z^ + Xi{l^n)] {l^a)z + {l + a)^/z^ + \i 

The term in the second bracket is non- vanishing for Re(z) > if and only if |a| < 1 (see Lemma 1 of Ref. |30] for a 
proof). The term in the first bracket never vanishes provided that O < 1. 

Summarizing, the CPBC with Weyl control and a choice of parameters that allows for fi = 1 or such that ( = — Ai 
and f2 < 1 does always yield an initial-boundary value formulation that satisfies the determinant condition, as long 
as the main evolution system is strongly hyperbolic and the constraint propagation system is symmetric hyperbolic 
(see sections Ull and 1111}) . 



B. CPBC without the Weyl control 

Next, we analyze the family of boundary conditions (|37I38|) where we set the boundary data h and Hab to zero. In 
this case, the result is less robust and depends more strongly on the choice of the parameters. For example, assume 
again that f2 = 1 which implies that — 0. Now Eq. I|69() has to be replaced by 



(1 + b){2z^ - C)Vz^ + l + {l - b)z{2z^ + 1) 



03 =0. 



By taking z real and positive and by taking the limits z — > and z — > cx3 one sees that the expression inside the 
square bracket changes its sign if — 1 < 6 < 1 and C > 0. In this case it does not follow that (T3 is zero and the system 
suffers from the presence of ill posed modes. 

For simplicity, let us assume that Ai = A2 = A3 = 1 and that ^ = — 1. All of the simulations below will satisfy these 
conditions. In those cases, we obtain 



/b(z)Vz2 + lai = 0, (70) 
fb{z) [(2z2 + 1)^3 + as] - 0, (71) 
fa{z) [(1 - r!)(a3 + ae) + z^og] = 0, (72) 



where fa(z) = (1 — a)z -I- (1 -t- a)\/ z^ + 1. Since \a\ < 1, \b\ < 1, fa and ft never vanish for Re(z) > and the only 
condition we obtain is < 3/2. 



VI. NUMERICAL IMPLEMENTATION 



In this section we discuss how to numerically implement the IB VP with or without Weyl control. For quasilinear first 
order symmetric hyperbolic systems with maximally dissipative boundary conditions there are well known methods 
for discretizing the problem such that numerical stability is guaranteed at the linearized level (see [sol Is?. 'EE. "sg] for a 
recent application in the context of numerical relativity, and references therein for the original work) . Unfortunately, 
in our case, the boundary conditions are more complicated since the CPBC (|35|) and the conditions 140|) which control 
the Weyl tensor are not even algebraic, and so we cannot apply this methods for discretizing the boundary conditions. 
So we first describe our method for implementing the constraint-preserving boundary conditions, and then explain 
how the differential equations are discretized. 

For simplicity, we focus on the case where the coupling coefficients Sj, a, b and c vanish, which corresponds to 
Sommerfeld-like conditions, as discussed in Section Hvl The boundary conditions (|35|l and H40() are implemented in 
the following way: For all points that lie on the boundary we add to the right-hand side of the evolution equations 
terms that are linear in the quantities V^^'^ and w^^b ~ Hab (which are zero if the boundary conditions are satisfied). 
More precisely, we replace the right-hand side of equations 13I4I5|) by 

doK,, = (as before) +p,°,.yW+g,^^(z«i+J-i?^s), (73) 
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dodktj = {as hefoie) + pl,jV^+^ 
doA, = (as before) +p°V;(+) + 



Ikij y^AB 



gf^(^W-i7^B), 



Hab), 



(74) 
(75) 



where the matrix coefficients p°-, Qij , etc. are allowed to depend on the three metric and the unit outward normal 
one-form to the boundary. (The q's are set to zero for the case without Weyl control.) It is clear that these extra 
terms change the principal part of the equations. The idea is to choose the p's and q's in such a way that with respect 
to the unit outward normal rii to the boundary the ingoing characteristic fields , and become zero speed 
fields while the speed of the other fields remains unchanged. This has the effect of eliminating the normal derivatives 
in the evolution equations for v^^J , and , and hence these equations are intrinsic to the boundary. We will 
see that this requirement uniquely determines the p's and q's. In order to see this we first notice that the ingoing 
characteristic constraint fields have the form 



V. 



(+) 



w 



(+) 
ab 



_l_d_ 



+ QAB{d\\u,u), 



Q{d\\u,u), 
+ QA{d\\u,u), 



(76) 
(77) 
(78) 



where d/dn denotes the normal derivative to the boundary and where the expressions Q, Qa and Qab only depend 
on the variables u and their derivatives d\\u tangential to the boundary. Therefore, the principal symbol corresponding 
to the modified system (|73I74I75|) and the direction of the unit outward normal Ui is given by 





= ±\ 


/^i^iV + 


(±)ay(+) , J±)CDA±) 
Pnn ''a ' Hnn '-'CD^ 


(±) 


= ±\ 


A^^Si + 


{±)ay(+) {±)CD (±) 
PaA +1aA ^CD^ 


(±) 


= ±\ 


/Aie; + 


(±)ayi+) i±)CD.{±) 
PnA +9nA ' 


-(±) 

^^^AB 


= ±V 


AB + PaB 


+1aB ^CD^ 


(0) 
f^^Ann 


= 04 


l^Ann ^ a 


) , JO)CD {±) 



(79) 

(80) 

(81) 

(82) 

(83) 
(84) 

where p^t\ ... PaIu and git\ ... (7^°^„ are defined in terms of iptj.plipPi ) and {qt, , qlij , qt) ^ ^'^^ same way as v'i^ , 



,(0) 



in terms of {Kij, dkij,Ai), and where 



^ 2 



(85) 
(86) 

Assume first that A2 = A3 = 1. In this case, we see immediately that win'', ... ^^Ann '"^^ only remain characteristic 
fields (i.e. fields with respect to which the principal symbol is diagonal) if all p's and q's are zero except for p^", 
p'nA^ and (fAB^^ ■ Then, the choice 



Paa 

PnA ^B 
-( + )CI5„.( + ) 
^AB 



W, 



CD 



(87) 



(+) 

AB 



(89) 



^AA' ^'nA and "C^. The matrix coefficients ... q^^ are easily obtained from 
this by applying the inverse transformation to the one that defines the characteristic fields in terms of the main 



yields zero speeds for the variables wVV, v'l^J and v 



variables (see |B)| . If A2 7^ 1 or A3 7^ 1 it is not possible to retain the fields 

„(+) 

case, we replace v^^l, v^^} by v^^^l = X2Kj^ + and v^^} = X^iKnA + DnA 
given in0 and set 



V. 



(+) 



as characteristic fields. In this 

,A 



AA-> 

respectively, where D\ and DnA are 



Paa 



VA2 



^ n ; 



(90) 
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^r^^^i = -A'l (92) 

and all other p's and g's to zero. This implies ^iv^^l — M^Ia'' = 0' ^^''^^AB — so and are zero speed 

fields while the speeds of the remaining characteristic fields are unchanged. 

Next, we discretize the evolution system (I1I2I73I74I75|I using the method of lines. The matrices p and q are set to 
zero at interior points and are chosen as described above at boundary points. The spatial discretization uses difference 
operators satisfying summation by parts (SBP) (see, for example, (431] ) . In this paper we use two of these difference 
aerators, which we call D2-1 and D8-4. These operators satisfy SBP with respect to diagonal norms; it can be seen 
that the use of these kind of norms implies that the order of the operator at and near boundary points is half 
that one in the interior. 

D2-1 is a simple operator which is second order accurate in the interior and first order accurate at boundaries. If 
the gridpoints range from to TV and the gridspacing is denoted by Aa;, this operator is given by 

Duq = -r^- (ui - Mo) , 
A.X 

Dui = (ui+i - Wi_i) for i = 1, . . . ,iV - 1, 

2Ax 

Dun = -J— {uN ~ UN-i) ■ 
Ax 

We typically add Kreiss-Oliger like dissipation |^ modified at and near boundaries such that the resulting operator is 
negative definite with respect to the norm for which SBP holds (much in the way Kreiss-Oliger dissipation is negative 
definite in the absence of boundaries) . This operator, which we call Qd, is added to the right-hand side of the evolution 
equations with a free (nonnegative) multiplicative parameter a diss (see [sSj for more details): 

QdUo = -2adiss^xD\u^, 

QdUi = -(Jdiss^x{D\ - 2D+D^)ui, 

QdUi = -cJd^ss{^xf{D+D-fui, fori = 2,...,iV-2, (93) 

QdUN-i = -(TdissAx^D'i - 2D+D-)uN-i, 

QdUN = -2(Tdiss^xD'^_UN- 



D8-4 is one of Strand's three-parametric difference operators of order eight in the interior and four at and close to 
boundaries, chosen in Ref. [60| to minimize its spectral radius and therefore its associated Courant limit. We have 
not been able to extend the associated Kreiss-Oliger dissipation for this operator in a straightforward way such as to 
include the presence of boundaries (see l60'| for more details). Therefore, here, we simply set the dissipative operator 
to zero near boundaries. However, as we will discuss later, the resulting dissipation operator does not seem to yield 
a stable numerical scheme in the presence of CPBC (as opposed to maximally dissipative boundary conditions) and 
a better dissipation operator is needed. 

For the time discretization we use a third order Runge-Kutta integrator. The boundary conditions (|37|) and (for 
the case without Weyl control) are imposed in the following way: We first compute the vector field corresponding 
to the discretized right-hand side of the system (|1I2I73I74I75|I and represent it in the characteristic basis. Then, 

the fields corresponding to Vm are overwritten by the time derivative of h and (for the case without Weyl control) 
the field corresponding to are overwritten with the time derivative of Hab- Finally, the resulting vector field 
is transformed back to the original basis. In all the simulations presented below we use the same gridspacing in all 
directions. Ax = Ay = Az and a Courant factor A of either A = 0.25 or A = 0.5. The code we use is based on the 
one described in Ref. [sof. 



VII. SIMULATIONS: SPACETIMES EVOLVED AND FORMULATIONS OF THE IB VP USED 



A. Spacetimes evolved 



The following initial data is used in the 3D evolutions of the next section. In each case, the spatial domain is the 
cube {(x, y, z) € [—1, 1]"^} of side length 2. 
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1. Random data or robust stability test 



Here we consider initial data corresponding to flat space and add a random perturbation to it. Therefore, initially, 
the fields are chosen to be 



9^, = % + (94) 

= en% (95) 

4,:, = (96) 

N = l + e7^" (97) 

A, = eTll (98) 



where the different TZ quantities are random numbers which are uniformly distributed in [—1, 1]. Similarly, at each 
timestep the boundary data h, Hab and Hab (see Eqs. (jSIJ, and respectively) for the "gauge" and 

"physical" degrees of freedom is set to a random number of the same order, eTZ. In the simulations shown below we 
choose e = 10~^. 

This robust stability test [45l |46| is designed to test the numerical stability of the scheme by finding out whether 
there is, at fixed time, unbounded growth as resolution is increased. One of the advantages of the test is that it 
excites all frequency modes allowed by a given resolution and it is therefore useful in spotting instabilities (if present) 
that could otherwise remain hidden in some convergence tests. However, this test is not a convergence test, since the 
constraints are not satisfied and since different random data is used when resolution is changed. 



2. Brill waves 



Axisymmetric Brill waves nWi are evolved in the next section. The corresponding initial data for the three-metric 
is given in Cartesian Coordinates x, y, z in the form 



ds^ = gijdx^dx^ = '^^ 



2q ( xdx + ydy 



+ e^'^dz^ + 



xdy — ydx 
P 



(99) 



where p ~ \J x^ + y^, the function q has the form 



q = Ap^ cxp — 



(100) 



and the conformal factor is obtained by solving the Hamiltonian constraint for time-symmetric initial data. In 
order to do so, we use the numerical elliptic solver BAM and the IDBrill Thorn, both publicly available from the 
CACTUS distribution 54]. The initial data for Kij, N, Ai and dkij is given by 



N = 
A^ = 

^kij 



0, 

1, 

0, 



(101) 
(102) 
(103) 
(104) 



where Dk is the D2-1 finite difference operator in the k direction. For these Brill wave simulations we choose the D2-1 
operator not only for computing the initial data for dkij but also for the discretization of the right-hand side of the 
evolution equations. There is no advantage in using a higher order accurate difference operator because the elliptic 
solver is only second order accurate. In the evolutions shown below, (T^ = 1/5. For weak enough waves the solution 
bounces at the origin and disperses to infinity, while for strong waves an apparent horizon is typically found |55|. In 
the evolutions below we use very weak waves, corresponding to an amplitude A = 10^^ (while the critical solution 
is believed to correspond to A ^ 4.8 [E^). As we shall see, even in these very weak field evolutions there is a clear 
violations of the constraints when non constraint-preserving boundary conditions are used. Even though these waves 
are axisymmetric, we have evolved them in full 3D. 
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3. Gauge solutions 

Initial data for a static solution corresponding to a gauge transformation of flat spacetime can be obtained by 
setting Kij, N and Ai to the same expressions as above [Eas. (|101llU2llUt^ ]. while the three-metric is obtained by 
starting from the flat metric in spherical coordinates, 

ds^ = dr^ + r^{d^^ + sm^i}dip^), 

performing a coordinates transformation 

r = f(l-ae-*^'/'^o) 

and transforming the resulting metric to the Cartesian coordinates (x,y,z) — (r sin "i? cos 0, f sin "i? sine/), f cos ■(?). Initial 
data for d^ij is obtained by computing analytically the gradient of the resulting Cartesian components of the three 
metric, d^ij = dugij- This gauge solution has been used to compare the stability properties of different formulations 
of Einstein's vacuum equations In our simulations, we choose a = 0.1 and ctq = 0.2. 

4. Injecting pulses of gravitational radiation through the boundaries 

Finally, we consider an example where we start with flat initial data, i.e. iV = 1, gij = 5ij^ Kij — 0, dkij — 0, 
Ak — 0, but inject a pulse of gravitational radiation by choosing as boundary data 

/i = 0, (105) 
h22 = = A23 = -2a(i - <o)^Tr'e[~^*-*")'/"' (106) 

where p is a "polar coordinate at each face". For example, at the a: = ±1 faces p — \/]p + z^ and similarly for the 
other faces. In the simulations shown below we choose = 1-2, (7f = cr^ = 0.2 and amplitudes a = 0.01 and 2. 

B. Formulations of the IB VP used 

In the simulations below, we choose F — NK corresponding to time harmonic slicing, set the shift to zero everywhere 
and at all times and restrict ourselves to the case in which the parameter Q in Eq. ^ is set to C = — 1. Furthermore, we 
only consider choices of parameters for which A2 = 1 and A3 = 1, and set the parameter deflned in Eq. H62|l to zero. 
This implies that the characteristic directions lie along the light cone or the normal to the t — const hypersurfaces. 

There are two subsets of parameter space which fulfill these requirements: 

1. Mono-parametric family (parametrized by the nonzero value of %): 

7 = -^, C = -l, r;-2, ^ = -f, Xy^O. (107) 

2. Bi-parametric family (parametrized by the parameters 77 and 7 7^ —1/2): 

C = -l, X^-^, ^ = -f+^-2> 7^4, (108) 

One can show that for these families the evolution system H1I2I3I4I5|I is symmetric hyperbolic. However, notice that 
the mono-parametric family violates the condition H28|) which means that the constraint propagation system is not 
symmetric hyperbolic. The bi-parametric family satisfies the conditions (|28I29I30() if and only if < 77 < 2. According 
to the analysis in the previous section the determinant condition is satisfied for this subfamily. In this article, we 
consider the following four cases: 

• CPBC without Weyl control: a completely ill posed case 
The bi-parametric family is used in this case, with parameters (77,7) = (3, 0), and CPBC without Weyl control. 
As mentioned above, the resulting parameters violate the condition for the constraint propagation system to be 
symmetric hyperbolic. Therefore, there is no guarantee that the IBVP is well posed. As a matter of fact, it 
turns out that the determinant condition is violated in this case and thus the problem admits ill posed modes. 
All the simulations performed below confirm this fact. 
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• CBPC with Weyl control 

Here we choose the bi-parametric family and CPBC with Weyl control. The parameters 77 and 7 are chosen 
to be either (ry,7) = (1,0) or (77,7) = (7/4,-2/3); both choices satisfy the determinant condition. Although 
the growth rate is smaller for the second choice, in both cases we find that the system is numerically unstable. 
However, we also notice that the instability seems to be milder than in the completely ill posed case: In the 
weak Brill wave runs, for instance, one needs to run for several crossing times or use very high resolution before 
noticing the lack of convergence. A similar situation arises in the numerical evolution of weakly hyperbolic 
systems 47] with periodic boundary conditions. 



• CPBC without Weyl control 

Again, we choose the bi-parametric family, but now we consider CPBC without Weyl control, and choose 
(77,7) = (1,0). The numerical simulations presented below suggest that the system is numerically stable. 

• Maximally dissipative 

The mono-parametric family with x = ^1 is used with maximally dissipative boundary conditions. This 
system should be stable, since the evolution equations are symmetric hyperbolic and so the IBVP is well 
posed. The simulations below confirm this. But more importantly, they provide an explicit demonstration 
that evolutions of a system with constraint-preserving boundary conditions, when numerically stable, are more 
accurate than standard, maximally dissipative ones, as in the latter case the boundary conditions introduce 
constraint violations that do not converge to zero. 



For each of the four IBVP described above we run the robust stability test and evolve weak Brill waves. For the 
case of CPBC without Weyl control, we also evolve the gauge solution and inject pulses of gravitational radiation 
from the boundary of an initially flat spacetime. 



For these runs, the Courant factor is A = 0.5 and the resolution varies from 21^^, 41'^ to 81"^ gridpoints. The spatial 
derivatives are discretized using the D2-1 operator, and a dissipation parameter of adiss — 0.05 is chosen. Recall 
that random initial and boundary data is given here. In Figure ^ we show the energy of the constraints versus light 
crossing time. The energy of the constraints is defined to be the sum over all gridpoints of the sum of the square of 
the components of the constraint variables divided by the number of gridpoints. 

As expected, in the completely ill posed case, the energy of the main variables grows at fixed time as resolution is 
increased, which makes the code crash in the timescale of less than a crossing time (there is also an increase in the 
energy at fixed resolution as a function of time, though this does not necessarily represent a numerical instability). 
Figure^shows that there is also a similar kind of growth in the energy for the constraints. The timescale and explosive 
kind of the numerical instability are similar to those found in initial value problems that are completely ill posed due 
to the the presence of complex eigenvalues in the principal part 5^ l6l|. This similarity is, indeed, expected, as at 
the analytical level the problem admits exponentially in time growing modes, where the exponential factor increases 
with the frequency as shown in section Recall, however, that the main evolution system is symmetric hyperbolic; 
thus the boundary conditions are responsible for the instabilities. 

In the CPBC case with Weyl control the runs also show evidence of a numerical instability both in the main 
and constraint variables, though the growth rate is somehow smaller than before. What is somehow unexpected is 
the resolution dependent growth in the constraint variables. After all, at the continuum, the constraints' evolution 
is governed by a symmetric hyperbolic system with maximally dissipative boundary conditions which constitutes a 
well posed problem by itself. However, it seems that the main system is unstable and that at the discrete level, 
the instabilities do have an effect on the constraint variables. In order to gain a better insight into this problem 
we have analyzed the semi-discrete problem, where only the space derivatives are discretized. Even in the simpler 
case of linearizations about Minkowski spacetime we found that while the discrete constraints do obey a symmetric 
hyperbolic system it is far from obvious that the way in which we implement the constraint-preserving boundary 
conditions allow for a semi-discrete energy estimate. It also seems difficult to represent the boundary conditions in a 
different way such that a semi-discrete energy estimate can be shown for the constraint propagation system. We will 
not attempt to address this question further in this article. Notice that even if we found such a discretization, the 
system could still suffer from the presence of more weakly ill posed gauge or physical modes that are undetected by 
the determinant condition. An explicit example of such weakly ill posed modes is given inIO 



VIII. SIMULATIONS 



A. Random, or robust stability test 
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FIG. 1: Random data runs for the D2-1 derivative, with Courant factor A — 0.5 and dissipation Cdiss = 0.05. The energy of the 
constraints is defined to be the sum over aU gridpoints of the sum of the square of the components of the constraint variables 
divided by the number of gridpoints. 



As opposed to the previous cases, in the CPBC case without Weyl control the runs strongly suggest that the 
resulting systems is numerically stable and that the continuum problem might be well posed. At early times the 
constraints' energy decrease while at around 100 crossing times the constraint variables start to slowly grow in time. 
However, the growth rate does not become larger with increasing resolution, as opposed to the previous cases. 

Finally, the random runs with maximally dissipative boundary conditions also strongly suggest that the system is 
numerically stable, though this is expected because in this case we know that the IBVP is well posed. 
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B. Brill wave evolutions 



In these simulations we use the smaller Courant factor of A = 0.25. Figure [21 shows weak Brill wave runs for the 
four formulations of the IB VP described in Section IVII Bl The expectations based on the random data runs of the 
previous section are confirmed. Namely, the completely ill posed CPBC case is manifestly numerically unstable. The 
CPBC with Weyl control case that satisfies the determinant condition seems to be unstable as well, though in a 
"weaker" sense while the two remaining cases (the maximally dissipative one and the CPBC without Weyl control 
case that satisfies the determinant condition) seem to be numerically stable. 
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FIG. 2: Weak Brill runs for the four formulations of the IBVP used in the simulations of Fig0 The expectations based on 
that figure regarding numerical stability (or its lack) are here confirmed. 



Figure |31 shows the last two cases on a shorter timescale and with a higher resolution run added. Notice how in the 
CPBC case the constraint variables seem to converge to zero, while in the maximally dissipative case with the same 
initial data and resolution the lack of convergence to zero is evident. 



C. Gauge solutions 

Here, we concentrate on the CPBC without Weyl control case. Since the initial data is given in analytic form we 
use the high order accurate finite differencing operator D8-4. The results are shown in Figure 01 Notice that the high 
order accurate scheme results in a much faster convergence of the constraints to zero as resolution is increased. Here, 
the Courant factor is A = 0.25 and the dissipation parameter a^iss = lO"**. 
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FIG. 3: Same runs as previous figure for tlie two numerically stable cases, but on a shorter timescale and with more resolution. 
Notice the lack of convergence to zero in the maximally dissipative case, as opposed to the constraint-preserving boundary 
conditions case. 
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FIG. 4: Gauge solution simulations with a high order scheme, Courant factor A = 0.25, and dissipation adiss = 10 



D. Injecting pulses of gravitational radiation through the boundaries 

As in the previous gauge simulations, we use the CPBC formulation without Weyl control. However, we have found 
that in this case for the same resolutions used in the previous simulations a numerical instability shows up after some 
time. This is probably due to the fact that the dissipative operator is zero near the boundary and is not negative 
semi-definite with respect to the scalar product for which SBP holds. Therefore for this case we have used the D2-1 
operator, with Courant factor A = 0.25, and dissipation parameter adiss = 0.05. In Figure [3 we show the energy of 
the constraints as a function of time for two resolutions and the two amplitudes a = 0.01 and a = 2. The energy 
starts from zero, as the initial data consists of Minkowski spacetime. After a short time the energy quickly grows as 
the pulses arc injected into the domain through the six boundaries, until the pulses have been completely injected at 
roughly one crossing time, time at which the energy stays approximately constant (at fixed resolution). 

Figure |51 shows the maximum (in the computational domain) of the curvature invariant J := RabcdR"'^'''^ = 8{E'^ — 
B^). This curvature invariant starts at zero as well and remains small (compared to one) for amplitude a — 0.01 
while for amplitude a = 2 it increases to very large values as the pulses are injected. For example, at i = 1.2 we have 
J « 3 X 10^. To have an idea of how strong the curvature is, recall that for the Schwarzschild solution J — ASrn^ /r^, 
where m is the mass of the black hole and r is the area radial coordinate. Therefore, a curvature of J w 3 x 10^ 
would correspond to being inside a black hole of mass to = 1 at r « 0.7. This curvature is being caused solely by the 
injection of pulses through the boundaries, showing that the latter are able to handle very non-linear dynamics. 
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FIG. 5: Injecting pulses through the boundaries. Here,the Courant factor is A = 0.25 and the dissipation parameter is 
o'disa = 10~*. The energy starts from roundoff as the initial data corresponds to flat spacetime. The energy curves for the 
highest resolutions have been scaled by a factor of four in a way that overlapping with the curves for to the coarser resolution 
corresponds to second order convergence. 
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FIG. 6: Same simulations as those shown in the previous figure, showing here the maximum (in the computational domain) 
of the curvature invariant RatcdR'^'"^'^ ■ The values achieved by this invariant in the simulations with amplitude a = 2 are 
comparable to values found inside a Schwarzschild black hole. This curvature is produced solely by the injection of pulses 
through the boundaries. 



IX. CONCLUSIONS 



We derived a family of outer boundary conditions for first order hyperbolic formulations of Einstein's field equations 
with live gauges. These boundary conditions have the property of being constraint-preserving in the sense that they 
guarantee that any smooth solution to the evolution equations subject to these boundary conditions automatically 
satisfy the constraints if so initially. Furthermore, we have discussed different possibilities for using constraint- 
preserving boundary conditions in order to control the gauge and physical degrees of freedom at the boundary. One 
of these possibilities (which we called CPBC with Weyl control) is attractive from a physical point of view since it 
provides boundary data to the Weyl scalars ^'o and ^"4, which at null infinity represent the in- and outgoing radiation. 
Furthermore, these scalars are gauge-invariant for linearizations about the Kerr metric, an approximation that should 
be good in many simulations provided the boundaries can be moved sufficiently far from the strong field region. 
We also discussed simpler ways (which we called CPBC without Weyl control) of controlling the physical degrees of 
freedom; however, their physical interpretation is less clear. 

Next, we analyzed the stability of the resulting IBVP by analytical and numerical means. By considering high- 
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frequency perturbations we obtained the corresponding frozen coefficient problem which is hnear and can be analyzed 
using Fourier-Laplace transformations. In this way, one obtains a determinant condition which is a necessary con- 
dition for the well posedness of the problem. The satisfaction of this condition restricts the freedom in the choice 
of parameters in the formulation and the boundary conditions. In particular, we found that the violation of the 
determinant condition leads to the presence of ill posed constraint- violating or gauge modes. These modes arc ill 
posed in the sense that they grow exponentially in time with an exponential factor that is unbounded. One could 
say that the boundary conditions are responsible for the presence of these modes since, in the absence of boundaries, 
the initial value problem is well posed because the evolution equations are strongly hyperbolic. A further example of 
boundary conditions leading to an instability of an otherwise stable initial- value problem is given in [HI 

Next, we performed three-dimensional numerical simulations. In order to do so, we extended an earlier finite- 
differencing code and implemented the constraint-preserving boundary conditions. We first performed a robust sta- 
bility test which consists in specifying random initial and boundary data for different resolutions and checking that the 
time evolution of the fields does not exhibit unbounded resolution-dependent growth. We found that the set of CPBC 
without Weyl control considered in this article is robustly stable in this sense. We also compared the results from 
these boundary conditions with results obtained by simply freezing the ingoing fields of the main evolution system to 
their initial values. While both systems are robustly stable, we found very strong evidence for the freezing boundary 
condition to yield constraint-violating solutions: In contrast to the CPBC, the constraint variables seem to converge 
to a nonzero value for freezing boundary conditions. We further tested the CPBC without Weyl control by evolving a 
gauge solution and weak Brill waves, and found that in both cases the constraint variables seem to converge to zero. 
Finally, as an application, we considered a situation in which one starts with flat initial data and injects "pulses of 
gravitational radiation" through the boundaries, of enough amplitude to create very large curvature in the interior. 

We expect the CPBC without Weyl control to be useful for many applications, like improving the accuracy and 
stability of current binary black hole and binary neutron star simulations or for a successful implementation of 
characteristic or perturbative matching techniques 0j IS IS Q 

Unfortunately, our numerical results for the CPBC with Weyl control do not pass the robust stability test, although 
the instability in this case is much weaker than for the case of CPBC that violate the determinant condition. There 
are several possible causes for this instability. First, it is not clear if this instability is caused by a "bad" discretization 
of the problem - the continuum IB VP being well posed - or if the continuum problem is actually ill posed and thus the 
cause for the instability observed. Another possible source of problem is the cubical domain used in our simulations, 
which has edges and corners. It might be that we violate some compatibility conditions at such points and that this 
causes an instability [6^. Whatever the cause of the instability might be, it would be nice to have a better analytic 
insight into the problem. Although the determinant condition used here is necessary for the well posedness of the 
IB VP it is not sufficient, not even for the frozen coefficient problem. When making such a statement we have to 
specify in what sense we expect the IBVP to be well posed. In the linear regime, we expect the problem to be well 
posed with respect to a Hilbert space that controls the norm of the fields and the norm of the constraint 
variables. In this case, there are explicit examples [sS] that show that the satisfaction of the determinant condition 
is not sufficient for well posedness. A different examples which applies to the IBVP considered in this article is given 
in ini Clearly, it would be much more satisfactory to derive CPBC with and without Weyl control for which well 
posedness can be guaranteed. 
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APPENDIX A: A DERIVATION OF THE MAIN EVOLUTION AND CONSTRAINT PROPAGATION 

SYSTEM 

In this appendix we re-derive the main evolution system, carefully keeping track of the constraints that are being 
used, and find the constraint propagation system with the help of the (twice contracted) Bianchi identities. We start 
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with the 3 + 1 spht of the four-dimensional Ricci tensor, 

= (3)^^^. „ i.(3)y X3)v^.^ _ 2KaK'j + KK,j - d^K,, , (Al) 

where '^^'^Rij and '■^■'Vi are the Ricci tensor and the covariant derivative, respectively, belonging to the three metric 
Qij, N is the lapse, Kij the extrinsic curvature, and cJq denotes the operator N~^{dt — £f3)- In order to obtain first 
order equations we rewrite the spatial derivatives of the three-metric as 

dkgij = dkij — Ckij , (A2) 

where dkij are new fields and Ckij — are constraints. Correspondingly, we can split the Christoffel symbols, '^'^^F'^j^ , 
belonging to the three metric as 

(3)pfe _ -pfc _ -pfc /AQ^ 

where 

r% = i/' (2d(,,); - du,) , (A4) 

are Christoffel symbols belonging to a new connection V which is torsion-free but not necessarily metric (we have 
Vfc(7y = —Ckij, so V is only metric if the constraints Ckij — are satisfied), and where 



= -g'^' (2C(,,), - Cu,) . (A5) 



Since the symbols r\j correspond to the difference between the connections *^^^V and V, they are actually the 
components of tensor field. Substituting Eq. ljA2p into the expression for the Ricci tensor, we obtain 



= R,^ - R,^ , (A6) 

where Rij is given by Eq. H12(l and where 

Rij = ^k^% ^ ^{i^^j)k ~ ^\k^\j + ^\]^\k - ^\j'^\l + 2^(i'''^j)"= ~ 9''''Ck{ij}l ■ (A7) 

Here, we have defined Cikij = d[iCk]ij- Note that the hatted quantities vanish if Ckij = 0. Next, we replace the 
gradient of the logarithm of the lapse by a new field Ai minus a corresponding constraint variable. 



_ A MA) 



A-Cr, (A8) 

JV 

and rewrite 



2A,,df + cf^C^ + t%{Ak - Ci^'>). (A9) 



— - - vr ' vjiv ^ V (j/ij) -r i±i^j - V 

Using Eqs. (|X6|) and (|X9|) . we rewrite Eq. IjAip as 

doK.^j = i?y - V(,ylj) - A,Aj - 2K, 'isTy + KK^ - Ay , (AlO) 

where 

A,, ^ (4)^^^. + _ ^^cf - 2A(,c]f ) -I- Cf^cf ^ T\^{Ak - c[^^) (All) 

groups together the four Ricci tensor (which vanishes in vacuum) and the constraint variables. An evolution equation 
for Kij in vacuum can be obtained from the identity IjAlOp by setting A^j to zero. However, in order to obtain a 
strongly hyperbolic evolution system, one needs to set Ay equal to suitable combinations of the constraint variables 
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(see below). Next, in order to obtain evolution equations for the new fields dkij and Ai, we apply the operator do on 
Eqs. HA2IA8() and use the commutation relation 

+ ^ {Ts^,...^,.^k^,,P' + ... + T,,,,,„,^_,MJ') , (A12) 
for any r-rank symbol Ti-^i2,,A,., where the Lie derivative £fj of Ti-^i^_,i^ is formally defined by 

Using the evolution equations dogij — ~2Kij and d^N = —FlN, K, x'^) for the three metric and the lapse, respectively, 
we obtain 

2 

^odk^J = -2dkK,j - 2AkK,j + —g^^dj^dufi^ + Afe„ , (A13) 

where 

Afc,, = aoCfc,y + 2d^^K,, , (A15) 
A - f) r^A) dF >A) ^ 

Finally, we rewrite the Hamiltonian and momentum constraint. Let '^■'G^^ be the Einstein tensor, and let ^ — 
denote the contraction with the vector field do. Then, we have 

(4)Goo = C - C, (4)Go, = -(G, - G,), 

where the expressions for G and Cj are given by Eqs. H6I7|I and where 

G = ig^'i^M , Q = ^g^'iCki ~ 2CkuW , + \K'''Cjki ■ (A17) 

The main evolution equations are dogij = —2Kij, OqN = —F{N, K, x^), and Eqs. HA10IA13IA14|I where one sets the 
quantities 

Ay (7,0 = Ay +75«jG + C/'Gfc(,,)z , (A18) 
Afejj (?7,x) = J^ktj ~ V9k(iCj) - x9ijCk , (A19) 
A,(0 = A,-eG, , (A20) 

to zero. 

With this information it is not very difficult to find the constraint propagation system using the commutation 
relation ljA12|) and the twice contracted Bianchi identities (written in 3+1 form) 

ao^^^Goo = +^^'^ V [n'^^^Go.) + 2X(4)Goo + K^'^^^R^, - Kg^^^^^R,,, (A21) 
do^'^Goj = +^(3)V,(iV2(4)G„o)+if(^)Go 

+ {n^'^R.,) ~ ^(3)V, {Ng^^^^^Rrs) . (A22) 

Substituting Gqo = G — G, Goj = —{Cj — Cj) and using the equations Aij{-fX) = Oi ^kijiv^x) = and Ai{£_) = 
and Eqs. (|A11IA15IA16|I . a lengthy but straightforward calculation yields 

doC = - (l + X - I) g^'diCk + I.O., (A23) 

doCj = -(l + 27)ajG-5''5"a, (^Gfc[y,, + ^Gyfe,-CGfe(y)i^ 



^Oj 
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- rd,c[f + Z.O., (A24) 

9oCfey = Lo., (A25) 

9oC/Hj = 2 (.9i[/c5i]Ci +5i[fe5;]C0 +X5'y"5[,Cfe] +Lo., (A26) 

agC*^) = Lo., (A27) 

SoC^f^ = (A28) 

where we have defined C'l'^^ — N^^d[i{NCj^^) and where the lower order terms, l.o., are hnear algebraic expressions 
in the constraint variables, with coefRcients that depend on the main variables and their spatial derivatives. The 

(A) 

hard part of the calculation is to show that no spatial derivatives of Ckij and other than the antisymmetric ones 
(which can be re-expressed in terms of the constraint variables Cikij and C^^"^ , respectively) enter the lower order 
terms. 



APPENDIX B: CHARACTERISTIC FIELDS AND STRONG HYPERBOLICITY 



Given a first order evolution system with principal symbol A{n), the characteristic fields with respect to a fixed 
direction rii are defined to be the projections of the main fields onto the eigenspaces of A{n). In order to find the 
characteristic fields for the symbol defined by Eq. (|15|l . it is convenient to choose an orthonormal basis of three- vectors 
ei, 62, 63, such that e\ = = g^-'rij and express the main variables Kij, dkij and Ai with respect to this basis: 

Kab ~ ^ij^a^b ' ^cafc = ^kij^c^a^b ' ~ ^i^a ■ (^1) 

Here and in the following, we assume that is normalized such that g'^^niTij — 1. Assuming that Ai, A2 and A3 are 
positive (which is a necessary condition for strong hyperbolicity) and defining 

= — i — — ^— , if Ai 7^ At and Vl arbitrary otherwise, (B2) 

2(Ai - A2) ^ ' 

we can express the characteristic fields as {n refers to the the index a—\ and A, B, C, ... to the indices a = 2,3) 

I'lt^ = K„„ + nKi±^{Dr,„ + nDi), (B3) 



V 



(±) - i^A^ r,A 
AA 



Ki±^Di, (B4) 



V. 



(±) 

nA 



2 



= KnA ±^D„A, (B5) 



. A3 

(±) _ 



Kab ± Dab , (B6) 



where 



C 1 1 7 

Di = ]^5'^^[-dnAB + {l+0dABn]+l{hn-dn), 
DnA = -^(1 - OidnnA - ^a) + ^(1 + QdAnn - ^dA^ ^Aa , 
Dab = ^ [-dnAB + (1 + C)rf(AS)«] - J SabS^" [-dnCD + (1 + OdcDn] ■ 

Here, dk g^^dkij, bj = g'^'^dkij and Kab = Kab ~ SabS'~^^ Kcd/2. vim* , v^^J, v^^g have the speeds n = ±\/Ai, 
/i = ±\/A2, M = iV^S: a-iid fi = ±1, respectively. (These are the speeds with respect to the normal derivative 
operator. The speeds with respect to the time evolution vector field are obtained from these after the transformation 
fi t—^ Nfj, + P^Ui.) The remaining characteristic fields have speeds = 0. For the following representation, we assume 
that the conditions H28I29|I . which are necessary for the constraint propagation system to be symmetrizable, hold. 
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"Ann 



Defining w = (77 — 2% — 2) ^, the zero speed characteristic fields are 

N, (B7) 
9^j , (B8) 

dAnn - X^{bA - ^a), (B9) 
dAnB - ^ rjU}5AB{K - dn), (BIO) 

dABC - V^SA{B{bc) - dc)) - X^^^Bc{bA - dA), (Bll) 

A„ + acj [(2 + + 3x)6„ - (2r? + x)dn] - {2a + i)uj(bn - dn), (B12) 

AA~i^{hA-dA). (B13) 

Having obtained the characteristic fields explicitly, it is not difficult to verify the additional smoothness require- 
ments for the system to be strongly hyperbolic: Namely, we have to construct a symmetric positive definite matrix 
H = H(p, n, u) which symmetrizes the principal symbol A{p, n, u) defined in Eq. H15|) . The system is called strongly 
hyperbolic if H depends smoothly on p, n, u and is such that the matrix HA is symmetric for all p, n, u. The smooth- 
ness requirements is needed for the pseudo-differential calculus while the symmetry condition allows for appropriate 
energy estimates. The matrix H can be obtained from the quadratic form which is built by summing over the square 
of the eigenfields: 



(0) 
'^ABC 

„(0) 



Ann'^Bnn ' " " " AnC BnD 



Clearly, this quadratic form depends smoothly on rii since it can be written in terms of contractions with g^^ and n*. 
H is also smooth in the other variables provided that the parameters are smooth functions, and provided that the 
function Q can be chosen smoothly. The latter condition is nontrivial if there are crossing points in the values swept 
by Ai and A2. In all our simulations, we choose smooth parameters such that Ai = A2 = A3 = 1 and set = 0, so the 
system is strongly hyperbolic. 

Finally, we give the inverse transformation which allows to recover the main variables from the characteristic ones. 
First, compute 

(B14) 
(B15) 
(B16) 

- ^Ki , (B17) 



Kab = 


^[^AB- 


y ^AB 


KnA — 


2\:"nA " 


^ ^nA 


Ki ^ 


^(^AA 


y Vaa 


K — 


lfyi + ) _ 





from which Kab is obtained. Next, compute 

Dab 



DnA = (B19) 



Dl = ^{v'n-^\ll (B20) 



Dnn = ^{v'^^J - vU) ^Di . (B21) 



Next, set 



^" I \Dnn + Daa + i'i"' ; 



b \ ^ (2Daa-C5^ 
^" / \ Dnn + Daa 

"A \ — I Ann + " ^ABC (T>,')'^'\ 
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where 



/ 2+37-2Ljg+2Ljcr(27;+x-2) _ -, 



2(l+27+r)Cw) 
c^(»)+3x) + l 

1 



-1 



M2 = { , . (B25) 

Notice that as a consequence of the conditions (|28I29|) . 1 + 27 + r/^cj = — + 27)(2 + 2x — ^]) — rjQ < and the 
matrix Afi is weh defined. Using this, compute 

d-Ann = v^^l^+xi^{bA~dA), (B26) 

dAnB = VaIb + \ V^^ABibn - dn), (B27) 

dABC = Vabc + '^^^A{B{bc} - dc)) +Xi^SBcibA- dA), (B28) 

A„ = «i")-atj[(2 + r; + 3x)&n-(277 + x)rf«] + (2a + Ow(6n-rf„), (B29) 

= (B30) 

= -2bAB + {l + C)d(AB)n+SAB[7{bn-dn)-D^], (B31) 

^ ■ ' ■ ^ ■ ^dA, (B32) 



dnnA — C ^ 



2A.A - + C)'5 {dBCA - dABc) + ^ 

1 1 7 ' 

Dnn - - Qbn + - d« + A„ - - (6„ - d„) 



2 



(B33) 



from which one can compute the components of dkij and Ai with respect to the orthonormal basis. 

Finally, one obtains the coordinate components of the main variables by contracting with the dual basis 0°, which 
is defined by 0°(ej,) = (5;^ , a, 6 = 1, 2, 3: 

K,, ^ KabO'^^e) , dk^J = d^abOiet^'j , A, = AaO^ . (B34) 

APPENDIX C: A SPECIAL FAMILY OF SOLUTIONS 

In this appendix we show explicitly that the linearized IBVP with Weyl control cannot be well posed in if the 
parameter Vl defined in Eq. 1)621) is one; independent on whether or not the determinant condition is satisfied. In 
order to see this, we consider the following family of solutions: 

Let Wj be a one-form that satisfies d^Wj — 0, d^diWj — 0, and dtWj = 0, and let 

Kij = di^iWj] , (CI) 
dkij = -'^tdkd(.,Wj) , (C2) 
A, = 0. (C3) 

It is not difhcult to check that these expressions satisfy the evolution equations 144I45I46|I and the constraints d^{dk — 
bk) = 0, d[idk]ij = 0, d^Kij — djK — 0. For systems without boundaries, these solutions are trivial if appropriate 
fall off conditions on the fields are demanded since then the harmonic condition on Wj implies that it must vanish. 
However, if boundaries are present, Wj may be nontrivial. The electric and magnetic components of the linearized 
Weyl tensor corresponding to the solutions l|C2p are 



Eij — 0, Bij — -d(^iejysd'^w'^ , 

and Bij vanishes if the one- form Wj is closed. In particular, this is true if Wj is exact, i.e. if Wj = djf for some 
time-independent harmonic function /. In this case, we also have 



dlf ■ (C4) 
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Therefore, if = 1, the family of solutions ljC2p with Wj = djf and / harmonic and time-independent shows that the 
linearized IBVP with the boundary conditions (|35I37I40|I is not well posed in since then the boundary conditions 
are satisfied with homogeneous data and since the initial data depends only on second derivatives of / whereas dkij 
depends on third derivatives of / for i > 0. This results in frequency dependent growth of the solution of the form 
\k\t, where fc is a characteristic wave number of the initial data. 
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